

/** @file integrand.cc
 *
 * Implementation of integrand 
 *
 */

/*
 *  Copyright (C) 2007 Stefan Weinzierl
 *
 *  This program is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation; either version 2 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program; if not, write to the Free Software
 *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 */

#include <stdexcept>

#include "integrand.h"
#include "polyhedra.h"
#include "multi_index.h"
#include "constants.h"

namespace sector_decomposition {

  // -----------------------------------------------------------------------------------

  /**
   *
   * Constructor.
   *
   * The arguments are: arg_list_feynman_parameter is a vector containing the Feynman parameters as GiNaC::symbol, 
   * arg_epsilon is the expansion parameter, again of type GiNaC::symbol,
   * arg_order defines term in the Laurent series, which should be calculated.
   *
   */
  integration_data::integration_data(const std::vector<GiNaC::ex> & arg_list_feynman_parameter, GiNaC::ex arg_epsilon, int arg_order) :
    n(arg_list_feynman_parameter.size()), list_feynman_parameter(arg_list_feynman_parameter), 
    epsilon(arg_epsilon), order(arg_order)
  {}

  /**
   *
   * Output. Prints out the Feynman parameters, the expansion parameter and the order.
   *
   */
  std::ostream & operator<< (std::ostream & os, const integration_data & arg)
  {
    os << "class integration_data" << std::endl;

    os << " Feynman parameters: ";
    for (int j=0; j<arg.n; j++) os << arg.list_feynman_parameter[j] << " ";
    os << std::endl;

    os << " Expansion parameter: " << arg.epsilon << std::endl;

    os << " Order : " << arg.order << std::endl;

    return os;
  }

  // -----------------------------------------------------------------------------------

  /**
   *
   * Default constructor.
   *
   */
  exponent::exponent() : order_0(0), order_1(0)
  {}

  /**
   *
   * Constructor.
   *
   * order_0 is initialised with arg_order_0, 
   * order_1 is initialised with arg_order_1. 
   *
   */
  exponent::exponent(int arg_order_0, int arg_order_1) : order_0(arg_order_0), order_1(arg_order_1)
  {}

  /**
   *
   * Constructor.
   *
   * order_0 is initialised with arg_order_0, 
   * order_1 is initialised with arg_order_1. 
   *
   */
  exponent::exponent(GiNaC::ex arg_order_0, GiNaC::ex arg_order_1) : order_0(arg_order_0), order_1(arg_order_1)
  {
    if ( !GiNaC::is_a<GiNaC::numeric>(arg_order_0) ) throw std::runtime_error("Argument order_0 not numerical.");
    if ( !GiNaC::is_a<GiNaC::numeric>(arg_order_1) ) throw std::runtime_error("Argument order_1 not numerical.");

    if ( !(GiNaC::ex_to<GiNaC::numeric>(order_0).is_rational()) ) throw std::runtime_error("Argument order_0 not rational.");
    if ( !(GiNaC::ex_to<GiNaC::numeric>(order_1).is_rational()) ) throw std::runtime_error("Argument order_1 not rational.");
  }

  /**
   *
   * Subscription.
   *
   */
  GiNaC::ex exponent::operator[](size_t i) const
    {
      if (i==1) return order_1;

      return order_0;
    }

  /**
   *
   * Adds the argument to the current exponent.
   *
   */
  exponent & exponent::sum_up(const exponent & other)
    {
      order_0 += other.order_0;
      order_1 += other.order_1;

      return *this;
    }

  /**
   *
   * Subtracts the argument from the current exponent.
   *
   */
  exponent & exponent::subtract_off(const exponent & other)
    {
      order_0 -= other.order_0;
      order_1 -= other.order_1;

      return *this;
    }

  /**
   *
   * Adds lambda times the argument to the current exponent.
   *
   */
  exponent & exponent::sum_up(const exponent & other, int lambda)
    {
      order_0 += other.order_0*lambda;
      order_1 += other.order_1*lambda;

      return *this;
    }

  /**
   *
   * Subtracts lambda times the argument from the current exponent.
   *
   */
  exponent & exponent::subtract_off(const exponent & other, int lambda)
    {
      order_0 -= other.order_0*lambda;
      order_1 -= other.order_1*lambda;

      return *this;
    }

  /**
   *
   * Returns the integer part of the exponent.
   *
   */
  int exponent::get_integer_part() const
  {
    // returns the smallest integer larger or equal to order_0, if order_0 < 0
    // returns the largest integer smaller or equal to order_0, if order_0 > 0
    GiNaC::numeric a = GiNaC::ex_to<GiNaC::numeric>(order_0);
    GiNaC::numeric p = GiNaC::numer(a);
    GiNaC::numeric q = GiNaC::denom(a);

    GiNaC::numeric R = GiNaC::iquo(p,q);

    return R.to_int();
  }


  /**
   *
   * Returns the exponent as GiNaC::ex.
   *
   */
  GiNaC::ex exponent::get_exponent(const integration_data & global_data) const
  {
    return order_0 + order_1*global_data.epsilon;
  }

  /**
   *
   * Returns true if not zero.
   *
   */
  bool exponent::is_not_zero() const
  {
    if ( (order_0 != 0) || (order_1 != 0) ) return true;

    return false;
  }

  /**
   *
   * Output. Prints out the exponent in the form exponent(order_0, order_1).
   *
   */
  std::ostream & operator<< (std::ostream & os, const exponent & arg)
  {
    os << "exponent(" << arg.order_0 << "," << arg.order_1 << ")";

    return os;
  }

  // -----------------------------------------------------------------------------------

  /**
   *
   * Constructor.
   *
   */
  integrand::integrand() :
    index(), weight(), polynomial_list(), poly_exponent(), extra_factor(1)
  {}

  /**
   *
   * Constructor.
   *
   * index is initialised with arg_index,
   * weight is initialised with nu,
   * polynomial_list is initialised with poly_list,
   * poly_exponent is initialised with c.
   *
   * extra_factor is initialised with 1.
   *
   */
  integrand::integrand(const std::vector<size_t> & arg_index, const std::vector<exponent> & nu, 
		       const std::vector<GiNaC::ex> & poly_list, const std::vector<exponent> & c) :
    index(arg_index), weight(nu), polynomial_list(poly_list), poly_exponent(c), extra_factor(1)
  {}

  /**
   *
   * Constructor.
   *
   * weight is initialised with nu,
   * polynomial_list is initialised with poly_list,
   * poly_exponent is initialised with c.
   *
   * extra_factor is initialised with 1,
   * index is initialised with (0,1,2,...)
   *
   */
  integrand::integrand(const std::vector<exponent> & nu, 
		       const std::vector<GiNaC::ex> & poly_list, const std::vector<exponent> & c) :
    index(nu.size()), weight(nu), polynomial_list(poly_list), poly_exponent(c), extra_factor(1)
  {
    for (int k=0; k<nu.size(); k++) index[k]=k;
  }

  /**
   *
   * Constructor.
   *
   * polynomial_list is initialised with poly_list,
   * poly_exponent is initialised with c.
   *
   * extra_factor is initialised with 1,
   * index is initialised with (0,1,2,...,n-1),
   * weight is initialised with n copies of exponent(0,0).
   *
   */
  integrand::integrand(size_t n, 
		       const std::vector<GiNaC::ex> & poly_list, const std::vector<exponent> & c) :
    index(n), weight(n,exponent(0,0)), polynomial_list(poly_list), poly_exponent(c), extra_factor(1)
  {
    for (int k=0; k<n; k++) index[k]=k;
  }

  /**
   *
   * Constructor.
   *
   * index is initialised with arg_index,
   * weight is initialised with nu,
   * polynomial_list is initialised with poly_list,
   * poly_exponent is initialised with c.
   * extra_factor is initialised with arg_extra_factor.
   *
   */
  integrand::integrand(const std::vector<size_t> & arg_index, const std::vector<exponent> & nu, 
		       const std::vector<GiNaC::ex> & poly_list, const std::vector<exponent> & c,
		       GiNaC::ex arg_extra_factor) :
    index(arg_index), weight(nu), polynomial_list(poly_list), poly_exponent(c), extra_factor(arg_extra_factor)
  {}

  /**
   *
   * Constructor.
   *
   * index is initialised with arg_index,
   * weight is initialised with nu,
   * extra_factor is initialised with arg_extra_factor.
   *
   * polynomial_list and poly_exponent are empty.
   *
   */
  integrand::integrand(const std::vector<size_t> & arg_index, const std::vector<exponent> & nu, 
		       GiNaC::ex arg_extra_factor) :
    index(arg_index), weight(nu), polynomial_list(), poly_exponent(), extra_factor(arg_extra_factor)
  {}

  /**
   *
   * Returns the explicit form of the integrand.
   *
   */
  GiNaC::ex integrand::get_integrand(const integration_data & global_data) const
  {
    GiNaC::ex res = extra_factor;

    int n = index.size();
    for (int j=0; j<n; j++)
      {
	res *= pow( global_data.list_feynman_parameter[index[j]], weight[j].get_exponent(global_data) );
      }

    n = polynomial_list.size();
    for (int j=0; j<n; j++)
      {
	res *= pow( polynomial_list[j], poly_exponent[j].get_exponent(global_data) );
      }
  
    return res;
  }

  /**
   *
   * Returns the explicit form of the weight factor.
   *
   */
  GiNaC::ex integrand::get_weight_factor(const integration_data & global_data) const
  {
    GiNaC::ex res = 1;

    int n = index.size();
    for (int j=0; j<n; j++)
      {
	res *= pow( global_data.list_feynman_parameter[index[j]], weight[j].get_exponent(global_data) );
      }

    return res;
  }

  /**
   *
   * Returns the explicit form of the polynomials.
   *
   */
  GiNaC::ex integrand::get_polynomials(const integration_data & global_data) const
  {
    GiNaC::ex res = 1;

    int n = polynomial_list.size();
    for (int j=0; j<n; j++)
      {
	res *= pow( polynomial_list[j], poly_exponent[j].get_exponent(global_data) );
      }
  
    return res;
  }

  /**
   *
   * Returns the explicit form of the polynomials without the exponent.
   *
   * The polynomial is returned in expanded form.
   *
   */
  GiNaC::ex integrand::get_polynomials_without_exponent(const integration_data & global_data) const
  {
    GiNaC::ex res = 1;

    int n = polynomial_list.size();
    for (int j=0; j<n; j++)
      {
	res *= polynomial_list[j];
      }
  
    return res.expand();
  }

  /**
   *
   * Returns the prefactor.
   *
   */
  GiNaC::ex integrand::get_extra_factor() const
  {
    return extra_factor;
  }

  /**
   *
   * Returns true if the \f$k\f$'s polynomial is homogeneous.
   *
   */
  bool integrand::check_homogeneous_polynomial(size_t k, const integration_data & global_data) const
  {
    GiNaC::ex poly = polynomial_list[k];

    GiNaC::symbol x;

    size_t n = index.size();
    for (int j=0; j<n; j++)
      {
	poly = poly.subs( global_data.list_feynman_parameter[index[j]] == x );
      }

    if ( poly.degree(x) == poly.ldegree(x) ) return true;

    return false;
  }

  /**
   *
   * Returns the homogeneous degree, if the \f$k\f$'s polynomial is homogeneous.
   *
   * Returns \f$-1\f$ if it is not a homogeneous polynomial.
   *
   */
  int integrand::get_homogeneous_degree(size_t k, const integration_data & global_data) const
  {
    GiNaC::ex poly = polynomial_list[k];

    GiNaC::symbol x;

    size_t n = index.size();
    for (int j=0; j<n; j++)
      {
	poly = poly.subs( global_data.list_feynman_parameter[index[j]] == x );
      }

    int poly_degree = poly.degree(x);
    if ( poly_degree == poly.ldegree(x) ) return poly_degree;

    return -1;
  }

  /**
   *
   * Returns the degree of poly.
   *
   */
  int integrand::get_degree(GiNaC::ex poly, const integration_data & global_data) const
  {
    GiNaC::symbol x;

    size_t n = index.size();
    for (int j=0; j<n; j++)
      {
	poly = poly.subs( global_data.list_feynman_parameter[index[j]] == x );
      }

    return poly.degree(x);
  }

  /**
   *
   * Transforms the \f$k\f$'s polynomial into a homogeneous polynomial.
   *
   */
  void integrand::make_homogeneous(size_t k, const integration_data & global_data) 
  {
    GiNaC::ex poly = polynomial_list[k];
    int max_degree = get_degree(poly,global_data);

    GiNaC::ex res = 0;
    GiNaC::ex sum_parameters = 0;

    for (int j=0; j<global_data.n; j++) sum_parameters += global_data.list_feynman_parameter[j];

    poly = poly.expand();
    int n_ops = poly.nops();
    for (int j=0; j<n_ops; j++)
      {
	GiNaC::ex temp = poly.op(j);
	int local_degree = get_degree(temp,global_data);
	res += temp*pow(sum_parameters,max_degree-local_degree);
      }

    polynomial_list[k] = res;
  }

  /**
   *
   * Transforms all polynomials into homogeneous polynomials.
   *
   */
  void integrand::make_homogeneous(const integration_data & global_data) 
  {
    size_t r = polynomial_list.size();

    for (int k=0; k<r; k++)
      {
	if ( !check_homogeneous_polynomial(k, global_data) ) make_homogeneous(k, global_data);
      }
  }

  /**
   *
   * Restricts to the \f$k\f$'s primary sector.
   *
   */
  void integrand::restrict_to_primary_sector(size_t k, const integration_data & global_data)
  {
    size_t n = index.size();
    size_t r = polynomial_list.size();

    // get homogeneous degrees of the polynomials
    std::vector<int> hom_degree(r);
    for (int j=0; j<r; j++) hom_degree[j] = get_homogeneous_degree(j, global_data);
      
    // calculate degree
    exponent degree = exponent(-n,0);

    for (int j=0; j<n; j++) degree.subtract_off(weight[j]);

    for (int j=0; j<r; j++) degree.subtract_off(poly_exponent[j], hom_degree[j]);


    // substitute
    for (int j=0; j<r; j++)
      {
	polynomial_list[j] = (polynomial_list[j]).subs( global_data.list_feynman_parameter[index[k]] == 1 );
      }

    // remove variable x_k
    index.erase( index.begin()+k );

    weight.erase( weight.begin()+k );

    // add extra term
    if (degree.is_not_zero()) 
      {
	GiNaC::ex new_poly = 1;
	for (int j=0; j<index.size(); j++) new_poly += global_data.list_feynman_parameter[index[j]];

	polynomial_list.push_back(new_poly);
	poly_exponent.push_back(degree);
      }
  }

  /**
   *
   * Restricts to one primary sector for the integration over the standard hyper-cube.
   *
   */
  void integrand::restrict_to_primary_sector_cube(const multi_index_counter & index_set, const integration_data & global_data)
  {
    size_t n = index.size();
    size_t r = polynomial_list.size();

    // calculate degree
    exponent degree = exponent(-n,0);
    for (int j=0; j<n; j++) degree.subtract_off(weight[j]);
    extra_factor *= pow(2,degree.get_exponent(global_data));

    // loop over all Feynman parameters
    for (int k=0; k<n; k++)
      {
	GiNaC::ex x_k = global_data.list_feynman_parameter[index[k]];

	if (index_set[k]==0) // 0 < x_k < 1/2
	  {
	    // substitute
	    for (int j=0; j<r; j++)
	      {
		polynomial_list[j] = (polynomial_list[j]).subs( x_k == x_k/2 );
	      }
	  }
	else // index_set[k]==1,  1/2 < x_k < 1
	  {
	    // substitute
	    for (int j=0; j<r; j++)
	      {
		polynomial_list[j] = (polynomial_list[j]).subs( x_k == 1-x_k/2 );
	      }

	    // add polynomial
	    polynomial_list.push_back( 2-x_k );
	    poly_exponent.push_back( weight[k] );

	    // set exponent of x_k to zero
	    weight[k] = exponent(0,0);
	  }
      } // k=0; k<n;
  }

  /**
   *
   * Checks if the \f$k\f$'s polynomial vanishes if the elements of index_set are set to zero.
   *
   */
  bool integrand::check_zero(size_t k, const multi_index_ordered & index_set, const integration_data & global_data) const
  {
    GiNaC::ex poly = polynomial_list[k];

    size_t n = index_set.size();
    for (int j=0; j<n; j++)
      {
	poly = poly.subs( global_data.list_feynman_parameter[index[index_set[j]]] == 0 );
      }

    if ( poly == 0 ) return true;

    return false;
  }

  /**
   *
   * Returns true if at least one of the polynomials vanishes, when the elements of index_set are set to zero.
   *
   */
  bool integrand::check_zero(const multi_index_ordered & index_set, const integration_data & global_data) const
  {
    size_t r = polynomial_list.size();

    for (int j=0; j<r; j++) 
      {
	if ( check_zero(j, index_set, global_data) ) return true;
      }

    return false;
  }

  /**
   *
   * Returns the ldegree of the \f$k\f$'s polynomial in the elements of index_set.
   *
   */
  int integrand::get_ldegree(size_t k, const multi_index_ordered & index_set, const integration_data & global_data) const
  {
    GiNaC::ex poly = polynomial_list[k];

    GiNaC::symbol x;

    size_t n = index_set.size();
    for (int j=0; j<n; j++)
      {
	poly = poly.subs( global_data.list_feynman_parameter[index[index_set[j]]] == x );
      }

    return poly.ldegree(x);
  }


  /**
   *
   * Returns the ldegree of the product of all polynomials in the elements of index_set.
   *
   */
  int integrand::get_ldegree(const multi_index_ordered & index_set, const integration_data & global_data) const
  {
    size_t r = polynomial_list.size();

    int low_degree = 0;

    for (int j=0; j<r; j++) low_degree += get_ldegree(j, index_set, global_data);

    return low_degree;
  }

  /**
   *
   * Returns a polynomial, where only terms corresponding to corners of a polyhedra
   * have been kept.
   *
   */
  GiNaC::ex integrand::get_polyhedra_polynomial(const integration_data & global_data) const
  {
    GiNaC::ex poly = get_polynomials_without_exponent(global_data);

    // only one monomial ?
    if ( !GiNaC::is_a<GiNaC::add>(poly) ) return poly;

    size_t n = index.size();

    size_t j_poly = poly.nops();

    std::vector<bool> flag_corner(j_poly, bool(true));

    std::vector<GiNaC::ex> poly_vec;
    for (int j=0; j<j_poly; j++) poly_vec.push_back(poly.op(j));

    for (int j1=0; j1<j_poly; j1++)
      {
	GiNaC::ex monom_1 = poly_vec[j1];

	for (int j2=j1+1; j2<j_poly; j2++)
	  {
	    GiNaC::ex monom_2 = poly_vec[j2];

	    if ( flag_corner[j1] && flag_corner[j2] )
		{
		  int local_d_max = 0;
		  int local_d_min = 0;

		  for (size_t k=0; k<n; k++)
		  {
		    int d1 = monom_1.degree(global_data.list_feynman_parameter[index[k]]);
		    int d2 = monom_2.degree(global_data.list_feynman_parameter[index[k]]);

		    if (d1-d2 > local_d_max)
		      {
			local_d_max = d1-d2;
		      }
		    else if (d1-d2 < local_d_min)
		      {
			local_d_min = d1-d2;
		      }
		  } // k=0; k<n; k++
		
		  if ( local_d_min >= 0 )
		  {
		    flag_corner[j1] = false;
		  }
		  else if ( local_d_max <= 0 )
		  {
		    flag_corner[j2] = false;
		  }
		}
	  } // j2
      } // j1

    GiNaC::ex res = 0;
    for (int j=0; j<j_poly; j++) 
      {
	if ( flag_corner[j] ) res += poly_vec[j];
      }

    return res;
  }

  /**
   *
   * Implementation of Zeillinger's strategy.
   *
   */
  bool integrand::find_min_max(std::vector<size_t> & k_min_max_set, const integration_data & global_data)
  {
    bool flag_found = false;
    size_t k_min = 0;
    size_t k_max = 0;
    int length = 0;
    int count = 0;

    size_t n = index.size();

    GiNaC::ex poly = get_polyhedra_polynomial(global_data);

    if ( GiNaC::is_a<GiNaC::add>(poly) ) 
      {
	size_t j_poly = poly.nops();

	for (int j1=0; j1<j_poly; j1++)
	  {
	    GiNaC::ex monom_1 = poly.op(j1);

	    for (int j2=j1+1; j2<j_poly; j2++)
	      {
		GiNaC::ex monom_2 = poly.op(j2);

		int local_d_max = 0;
		int local_d_min = 0;
		size_t local_k_max, local_k_min;
		int local_max_count, local_min_count;

		for (size_t k=0; k<n; k++)
		  {
		    int d1 = monom_1.degree(global_data.list_feynman_parameter[index[k]]);
		    int d2 = monom_2.degree(global_data.list_feynman_parameter[index[k]]);

		    // there will be at least one component with d1-d2>0
		    if (d1-d2 > local_d_max)
		      {
			local_d_max = d1-d2;
			local_k_max = k;
			local_max_count = 1;
		      }
		    // and at least one component with d1-d2<0
		    else if (d1-d2 < local_d_min)
		      {
			local_d_min = d1-d2;
			local_k_min = k;
			local_min_count = 1;
		      }
		    else if ( d1-d2 == local_d_max )
		      {
			local_max_count++;
		      }
		    else if ( d1-d2 == local_d_min )
		      {
			local_min_count++;
		      }
		
		  } // k
		
		int local_length = local_d_max - local_d_min;
		int local_count = local_max_count + local_min_count;

		if ( (local_length < length) || ( (local_length == length) && (local_count < count) ) || (!flag_found) )
		  {
		    length = local_length;
		    count = local_count;
		    k_min = local_k_min;
		    k_max = local_k_max;
		    flag_found = true;
		  }

	      } // j2
	  } // j1

	if ( flag_found )
	  {
	    // std::cout << "length " << j_poly << " " << length << " " << count << std::endl;

	    if ( k_min < k_max )
	      {
		k_min_max_set.push_back(k_min);
		k_min_max_set.push_back(k_max);
	      }
	    else
	      {
		k_min_max_set.push_back(k_max);
		k_min_max_set.push_back(k_min);
	      }
	  }
      } // is_a<add>

    return flag_found;
  }

  /**
   *
   * Determines a subset of Feynman parameters accoriding to Spivakovsky's strategy.
   *
   */
  bool integrand::get_S_strategy_B(std::vector<size_t> & index_set, const integration_data & global_data)
  {
    GiNaC::ex poly = get_polynomials_without_exponent(global_data);

    bool flag = true;
    bool flag_degree_zero = false;

    std::vector<polyhedra> Delta;
    std::vector<polyhedra> Delta_tilde;

    std::vector<std::vector<size_t> > I; 

    std::vector<GiNaC::ex> order_string;

    I.push_back(index);

    Delta.push_back( polyhedra(poly,index,global_data).keep_corners_only() );
    if ( Delta.back().number_corners() == 0 ) flag = false;

    // iteration
    while ( flag )
      {
	// get Delta_tilde
	Delta_tilde.push_back( Delta.back() );
	Delta_tilde.back().factorize_single_variables();

	GiNaC::ex f = Delta_tilde.back().get_factorization_degree();

	order_string.push_back(f);

	if ( f == 0 ) 
	  {
	    flag_degree_zero= true;
	    flag = false;
	  }

	if ( flag )
	  {
	    // get S
	    std::vector<bool> S_bool;
	    Delta_tilde.back().find_S(S_bool, f);

	    // update I
	    std::vector<size_t> I_new;
	    for (size_t k=0; k<S_bool.size(); k++)
	      {
		if ( !(S_bool[k]) ) I_new.push_back( I.back()[k] );
	      }
	    I.push_back(I_new);
	    order_string.push_back( GiNaC::numeric(I_new.size()) );

	    // projection
	    polyhedra p1 = Delta_tilde.back();
	    p1.projection(S_bool,f,1);

	    polyhedra p2 = Delta.back();
	    p2.projection(S_bool,1,1);

	    Delta.push_back( polyhedra(p1,p2).keep_corners_only() );
	    if ( Delta.back().number_corners() == 0 ) 
	      {
		flag = false;
		order_string.push_back(Infinity);
	      }
	  }

      } // end iteration

    // last element of order_string
    order_string.push_back(Delta.back().get_factorization_degree() );

    if ( flag_degree_zero )
      {
	std::vector<size_t> Gamma;
	Delta.back().find_minimal_subsector(Gamma,1);

	for (size_t k=0; k<index.size(); k++)
	  {
	    bool flag_not_I_r= true;
	    for ( size_t j=0; j<I.back().size(); j++ )
	      {
		if ( index[k] == I.back()[j] ) flag_not_I_r = false;
	      }

	    bool flag_Gamma = false;
	    for ( size_t j=0; j<Gamma.size(); j++ )
	      {
		if ( index[k] == I.back()[Gamma[j]] ) flag_Gamma = true;
	      }

	    if ( flag_not_I_r || flag_Gamma ) index_set.push_back(k);
	  }
      }
    else // Delta_r = empty_set
      {
	for (size_t k=0; k<index.size(); k++)
	  {
	    bool flag_not_I_r= true;
	    for ( size_t j=0; j<I.back().size(); j++ )
	      {
		if ( index[k] == I.back()[j] ) flag_not_I_r = false;
	      }
	    if ( flag_not_I_r ) index_set.push_back(k);
	  }
      }

    // monomial ideal ?
    if ( flag_degree_zero && (Delta_tilde.size() == 1) ) return false;

    return true;
  }

  /**
   *
   * Determines a subset of Feynman parameters accoriding to strategy C.
   *
   */
  bool integrand::get_S_strategy_C(std::vector<size_t> & index_set, const integration_data & global_data)
  {
    GiNaC::ex poly = get_polynomials_without_exponent(global_data);

    bool flag = true;
    bool flag_degree_zero = false;

    std::vector<polyhedra> Delta;
    std::vector<polyhedra> Delta_tilde;

    std::vector<std::vector<size_t> > I; 

    std::vector<GiNaC::ex> control;

    std::vector<GiNaC::ex> order_string;

    control.push_back(0);
    I.push_back(index);

    Delta.push_back( polyhedra(poly,index,global_data).keep_corners_only() );
    if ( Delta.back().number_corners() == 0 ) flag = false;

    // iteration
    while ( flag )
      {
	// get Delta_tilde
	Delta_tilde.push_back( Delta.back() );
	Delta_tilde.back().factorize_single_variables();

	GiNaC::ex f_minus = control.back();
	GiNaC::ex f = Delta_tilde.back().get_factorization_degree();

	order_string.push_back(f);

	if ( f == 0 ) 
	  {
	    flag_degree_zero= true;
	    flag = false;
	  }

	if ( flag )
	  {
	    control.push_back(f);

	    // get S
	    std::vector<bool> S_bool;
	    Delta_tilde.back().find_S(S_bool, f);

	    int number_S_true = 0;
	    for (size_t k=0; k<S_bool.size(); k++)
	      {
		if ( S_bool[k] ) number_S_true++;
	      }

	    order_string.push_back( GiNaC::numeric( S_bool.size() - number_S_true ) );

	    // choose a subset of S_bool: here we always choose S_bool itself

	    // update I
	    std::vector<size_t> I_new;
	    for (size_t k=0; k<S_bool.size(); k++)
	      {
		if ( !(S_bool[k]) ) I_new.push_back( I.back()[k] );
	      }
	    I.push_back(I_new);

	    // projection
	    polyhedra p1 = Delta_tilde.back();
	    if ( f < f_minus )
	      {
		std::vector<GiNaC::ex> omega;
		Delta.back().get_omega(omega);
		for (size_t k=0; k<omega.size(); k++) omega[k] *= f/(f_minus-f);

		p1.add_point(omega);
	      }
	    p1.projection(S_bool,f,f);

	    Delta.push_back( p1.keep_corners_only() );

	    if ( Delta.back().number_corners() == 0 ) 
	      {
		flag = false;
		order_string.push_back(Infinity);
	      }
	  }

      } // end iteration

    // last element of order_string
    order_string.push_back(Delta.back().get_factorization_degree() );

    if ( flag_degree_zero )
      {
	std::vector<size_t> Gamma;
	Delta.back().find_minimal_subsector(Gamma, control.back() );

	for (size_t k=0; k<index.size(); k++)
	  {
	    bool flag_not_I_r= true;
	    for ( size_t j=0; j<I.back().size(); j++ )
	      {
		if ( index[k] == I.back()[j] ) flag_not_I_r = false;
	      }

	    bool flag_Gamma = false;
	    for ( size_t j=0; j<Gamma.size(); j++ )
	      {
		if ( index[k] == I.back()[Gamma[j]] ) flag_Gamma = true;
	      }

	    if ( flag_not_I_r || flag_Gamma ) index_set.push_back(k);
	  }
      }
    else // Delta_r = empty_set
      {
	for (size_t k=0; k<index.size(); k++)
	  {
	    bool flag_not_I_r= true;
	    for ( size_t j=0; j<I.back().size(); j++ )
	      {
		if ( index[k] == I.back()[j] ) flag_not_I_r = false;
	      }
	    if ( flag_not_I_r ) index_set.push_back(k);
	  }
      }

    // monomial ideal ?
    if ( flag_degree_zero && (Delta_tilde.size() == 1) ) return false;

    return true;
  }

  /**
   *
   * Creates the \f$k\f$'th subsector by remapping the variables given by index_set.
   *
   */
  void integrand::create_subsector(size_t k, const std::vector<size_t> & index_set, const integration_data & global_data)
  {
    size_t n = index.size();
    size_t r = polynomial_list.size();
    size_t m = index_set.size();


    // remap m-1 variables
    GiNaC::ex x_k = global_data.list_feynman_parameter[index[index_set[k]]];
    for (int j=0; j<m; j++)
      {
	if ( j != k )
	  {
	    GiNaC::ex x_j = global_data.list_feynman_parameter[index[index_set[j]]];

	    for (int l=0; l<r; l++)
	      {
		polynomial_list[l] = (polynomial_list[l]).subs( x_j == x_k*x_j );
	      }
	  }
      }

    // calculate degree of variable x_k
    exponent degree = exponent(m-1,0);

    for (int j=0; j<m; j++)
      {
	degree.sum_up( weight[index_set[j]] );
      }
    weight[index_set[k]] = degree;
  }

  /**
   *
   * Creates the \f$k\f$'th subsector by remapping the variables given by index_set.
   *
   */
  void integrand::create_subsector(size_t k, const multi_index_ordered & index_set, const integration_data & global_data)
  {
    size_t n = index.size();
    size_t r = polynomial_list.size();
    size_t m = index_set.size();


    // remap m-1 variables
    GiNaC::ex x_k = global_data.list_feynman_parameter[index[index_set[k]]];
    for (int j=0; j<m; j++)
      {
	if ( j != k )
	  {
	    GiNaC::ex x_j = global_data.list_feynman_parameter[index[index_set[j]]];

	    for (int l=0; l<r; l++)
	      {
		polynomial_list[l] = (polynomial_list[l]).subs( x_j == x_k*x_j );
	      }
	  }
      }

    // calculate degree of variable x_k
    exponent degree = exponent(m-1,0);

    for (int j=0; j<m; j++)
      {
	degree.sum_up( weight[index_set[j]] );
      }
    weight[index_set[k]] = degree;
  }

  /**
   *
   * The routine fills a vector, indicating how many subtractions are needed for each variable.
   * The indices for the corresponding variables are stored in index_set.
   * If no subtractions are needed, the corresponding entry is not included.
   * The return value is the number of variables, in which subtractions are needed.
   *
   */
  int integrand::get_subtraction_vector(std::vector<size_t> & subtr_vec, std::vector<size_t> & index_set) const
  {
    size_t n = index.size();

    for (int j=0; j<n; j++)
      {
	int integer_part = (weight[j]).get_integer_part();

	if ( integer_part < 0 ) 
	  {
	    index_set.push_back(j);
	    subtr_vec.push_back(-integer_part);
	  }
      }

    return index_set.size();
  }

  /**
   *
   * Returns one term in the Taylor expansion.
   *
   */
  GiNaC::ex integrand::taylor_term(GiNaC::ex expr, const multi_index_counter_indv & counter, const std::vector<size_t> & index_set, const integration_data & global_data) const
  {
    size_t m = index_set.size();

    GiNaC::ex res = expr;

    for (int j=0; j<m; j++)
      {
	GiNaC::ex x = global_data.list_feynman_parameter[index[index_set[j]]];
	res = res.diff(GiNaC::ex_to<GiNaC::symbol>(x), counter[j] );
	res = res.subs(x == 0);
	res = pow(x,counter[j])/GiNaC::factorial(counter[j]) * res;
      }

    // include sign
    return pow(GiNaC::numeric(-1),m)*res;
  }

  /**
   *
   * Returns one term in the Taylor expansion, integrated over the singular variables.
   *
   */
  void integrand::integrated_taylor_term(const multi_index_counter_indv & counter, const std::vector<size_t> & index_set, const integration_data & global_data) 
  {
    size_t m = index_set.size();

    GiNaC::ex res = extra_factor;

    for (int j=0; j<m; j++)
      {
	GiNaC::ex x = global_data.list_feynman_parameter[index[index_set[j]]];
	res = res.diff(GiNaC::ex_to<GiNaC::symbol>(x), counter[j] );
	res = res.subs(x == 0);
	res = res/GiNaC::factorial(counter[j])/( (weight[index_set[j]]).get_exponent(global_data) + counter[j] + 1 );
      }

    // include sign
    extra_factor = pow(GiNaC::numeric(-1),m+1)*res;

    // remove integrated variables
    for (int j=m-1; j>=0; j--)
      {
	index.erase( index.begin()+index_set[j] );
	weight.erase( weight.begin()+index_set[j] );
      }
  }

  /**
   *
   * Output. Prints out the data of the class integrand without reference to a class integration_data.
   *
   */
  std::ostream & operator<< (std::ostream & os, const integrand & arg)
  {
    os << "class integrand" << std::endl;

    os << " Index: ";
    for (int j=0; j<arg.index.size(); j++) os << arg.index[j] << " ";
    os << std::endl;

    os << " Weight: ";
    for (int j=0; j<arg.weight.size(); j++) os << arg.weight[j] << " ";
    os << std::endl;

    if ( arg.polynomial_list.size() > 0 ){
    os << " Polynomials: " << std::endl;
    for (int j=0; j<arg.polynomial_list.size(); j++) os << arg.polynomial_list[j] << ", " << arg.poly_exponent[j] << std::endl;
    }
    os << " Extra : " << arg.extra_factor << std::endl;

    return os;
  }

} // namespace sector_decomposition

